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We present a general method for obtaining the exact static solutions and collective excitation 
frequencies of a trapped Bose-Einstein condensate (BEC) with dipolar atomic interactions in the 
Thomas-Fermi regime. The method incorporates analytic expressions for the dipolar potential 
of an arbitrary polynomial density profile, thereby reducing the problem of handling non-local 
dipolar interactions to the solution of algebraic equations. We comprehensively map out the static 
solutions and excitation modes, including non-cylindrically symmetric traps, and also the case of 
negative scattering length where dipolar interactions stabilize an otherwise unstable condensate. 
The dynamical stability of the excitation modes gives insight into the onset of collapse of a dipolar 
BEC. We find that global collapse is consistently mediated by an anisotropic quadrupolar collective 
mode, although there are two trapping regimes in which the BEC is stable against quadrupolc 
fluctuations even as the ratio of the dipolar to s-wave interactions becomes infinite. Motivated by 
the possibility of fragmented BEC in a dipolar Bose gas due to the partially attractive interactions, 
we pay special attention to the scissors modes, which can provide a signature of superfluidity, and 
identify a long-range restoring force which is peculiar to dipolar systems. As part of the supporting 
material for this paper we provide the computer program used to make the calculations, including 
a graphical user interface. 

PACS numbers: 03.75.Kk, 34.20.Cf 



INTRODUCTION 



Since the realization of atomic Bose-Einstein conden- 
sates (BECs) in 1995 [I], there has been a surge of in- 
terest in quantum degenerate gases [21 [3] . Despite the 
dilutcness of these gases, interatomic interactions play 
an important role in determining their properties. In the 
majority of experiments the dominant interactions have 
been isotropic and asymptotically of the van der Waals 
type, falling off as 1/r 6 . At ultracold temperatures this 
leads to essentially pure s-wave scattering between the 
atoms. An exception to this rule is provided by gases 
that have significant dipole-dipole interactions [IHZI- In 
comparison to van der Waals type interactions, dipolar 
interactions are longer range and anisotropic, and this 
introduces rich new phenomena. For example, a series of 
experiments that have revealed the anisotropic nature of 
dipolar interactions are those on 52 Cr BECs in an exter- 
nal magnetic field. These have demonstrated anisotropic 
expansion of the condensate depending on the direction 
of polarization of the atomic dipoles [H |9] , collapse and 
rf-wave explosion |10j . and an enhanced stability against 
collapse in flattened geometries [IT]. Meanwhile, an ex- 
periment with 39 K atoms occupying different sites in a 
ID optical lattice has demonstrated the long-range na- 
ture of dipolar interactions in BECs through dephasing of 
Bloch oscillations [B] . Dipolar interactions have also been 
shown to be responsible for the formation of a spatially 
modulated structure of spin domains in a Rb spinor 
BEC 0. 



In order to incorporate atomic interactions into the 
Gross-Pitaevskii theory for the condensate one should use 
a pseudo-potential 2 , 3J . In the presence of both dipolar 
and van der Waals interactions the pseudo-potential can 
be written as the sum of two terms U(r) — £/ s (r) + ?7dd( r ) 
[T2lU5| . where r is the relative interatomic separation. 
The long-range dipolar interaction can be treated accu- 
rately within the Born approximation providing one is 
not close to a scattering resonance [21 [J3] . This first- 
order approximation means that the effective interaction 
is replaced by the potential itself. This is quite different 
to the shorter range van der Waals interaction, for which 
the Born approximation is not valid at low temperatures, 
and where one rather uses the contact potential 

U s (r) = g6(r). (1) 

The coupling constant g = 4irh 2 a s /m is given in terms 
of an s-wave scattering length a s and the atomic mass 
m. For the dipolar interaction, we consider two atoms 
whose dipoles are aligned by an external field pointing 
along the direction specified by the unit vector e. The 
potential is then given by 

^d(r) = %e i e,fc^l. (2) 



where C^d parameterizes the strength of the dipolar in- 
teractions, f is a unit vector in the direction of r, and 
summation over repeated indices is implied. A key figure 
of merit is the ratio of the two coupling strengths, defined 
as [Tfj] . 



£dd — Cdd/3.g. 



(3) 
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Dipole-dipole interactions can be either magnetic or elec- 
tric in origin. To date, the dipolar interactions seen in 
ultracold atom experiments [1H7] have all been magnetic 
dipolar interactions, for which Cdd = Mo^ 2 , where d is 
the magnetic dipole moment and fj,o is the permeabil- 
ity of free space. In terms of the Bohr magneton /is, the 
magnetic dipole moment of a 52 Cr atom is d — 6/ib giving 
£dd ~ 0.16 [4 . Although this is 36 times larger than the 
typical value of £dd found in the alkalis, it is still small. 
Thus, unless the system is in a configuration that makes 
it particularly sensitive [BJ, and/or is specially prepared 
[7] , the magnetic dipolar interactions in the atomic gases 
made so far tend to be masked by stronger s-wave inter- 
actions. In order to make dipolar interactions in BECs 
more visible, the Stuttgart group have succeeded in im- 
plementing magnetic Feshbach resonances [TT] in 52 Cr 
[TT] . These allow g to be tuned from positive to negative 
and even to zero. Moreover, the sign and amplitude of 
the effective value of Cdd can also be tuned by rapidly ro- 
tating the external polarizing field |16j ■ Polar molecules 
can have huge electric dipole moments and these systems 
are now close to reaching degeneracy [T5H2"3"] . By appro- 
priately tuning an external electric field a large degree 
of control can be exerted over these systems [23] . Com- 
bined with what has already been achieved in 52 Cr, one 
can realistically explore a large parameter space of inter- 
actions. 

The ground state of a trapped dipolar BEC has already 
been investigated theoretically by a number of authors, 
e.g. [T2HI51 l24Tl31| . with most studies focussing on the 
regime where g > and Cdd > 0. The presence of dipo- 
lar interactions was widely predicted to lead to certain 
distinctive effects, some of which have recently been seen 
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FIG. 1: Schematic illustration of the basic collective modes 
under consideration: the dipole mode D (shown here in the 
^-direction D x ), scissors mode Sc (shown here in x — z plane 
Sc xz ), the monopole mode M and the quadrupole modes Qi 
and Q2- These modes are discussed in more detail in Section 

eh 



experimentally. For example, if the dipoles are aligned 
in the ^-direction, then a condensate will elongate along 
z and become more "cigar" -shaped, i.e. undergo mag- 
netostriction, in order to benefit energetically from the 
attractive end-to-end interaction of dipoles. As £dd is 
increased, for example by reducing g with a Feshbach 
resonance, the BEC eventually becomes unstable to col- 
lapse, and this striking behavior has been realized in the 
experiment [TH]. Conversely, a condensate that is flat- 
tened by strong trapping along z will be mostly composed 
of repulsive side-by-side dipoles and so this "pancake" - 
shaped geometry is more stable, as confirmed experimen- 
tally [TJJ. In the limit that £dd becomes large, but the 
BEC remains in the pancake configuration due to tight 
trapping, remarkable density wave structures have been 
predicted for certain regions of parameter space close to 
the collapse threshold [28H30] . 

In this paper we work in the Thomas-Fermi (TF) 
regime, which is of rather general interest because it 
is formally equivalent to the hydrodynamic regime of 
zero-temperature superfluids [32] • The TF regime may 
be viewed as the semiclassical approximation to the full 
Gross-Pitaevskii theory. A stationary condensate enters 
the TF regime when the zero-point kinetic energy of the 
atoms due to the confinement by the trap becomes negli- 
gible in comparison to the total interaction and trapping 
energies. For BECs with repulsive interactions in a har- 
monic trap this generally occurs in the large N limit, 
where N is the number of atoms. However, for dipolar 
BECs the picture is considerably complicated by the par- 
tially attractive and partially repulsive nature of the in- 
teractions. The question of the validity of the TF regime 
in dipolar BECs has been addressed in [33] . 

The TF regime is theoretically simpler to handle than 
the full Gross-Pitaevskii theory, thereby facilitating an- 
alytical results. For example, under harmonic trapping 
it can be shown that the exact density profile of a dipo- 
lar condensate in the TF regime is an inverted parabola 
[231 |2T)] , similar to the usual s-wave case but distorted by 
the magnetostriction. Furthermore, the stability of the 
ground state to collapse can be estimated simply in the 
TF regime and reasonable agreement with experiment 
has been reported [TT] . Rotational instabilities of dipo- 
lar BECs are also amenable to analysis in the TF regime 
[3TI 135] . The current paper builds on these earlier works 
by applying the exact results available in the TF regime 
to collective excitations. 

The excited states of a BEC can be accurately calcu- 
lated within the TF regime provided they are of suffi- 
ciently long wavelength. The most basic collective exci- 
tations of a trapped BEC are the dipole (centre-of-mass) , 
monopole (breathing), quadrupole and scissors modes, il- 
lustrated schematically in Fig. [T] Their characterization 
offers important opportunities to measure interaction ef- 
fects, test theoretical models, and even detect weak forces 
|36j . Specifically, the scissors mode provides an impor- 
tant test for superfluidity [5TM40| . while the quadrupole 
mode plays a key role in the onset of vortex nucleation 
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in rotating condensates [3H [351 l4"TM45| . An instability 
of the quadrupole mode is also thought to be the mecha- 
nism by which collapse of dipolar BECs proceeds when it 
occurs globally Q31 [11 [3TJ S5] (rather than locally [3"T]). 
While the collective modes of a dipolar BEC have been 
studied previously [24j E3 EH I47H50] , key issues remain 
at large, for example, the regimes of Cdd < and g < 0, 
and the behaviour of the scissors modes. This provides 
the motivation for the current work. 

In this paper we present a general and accessible 
methodology for determining the static solutions and 
excitation frequencies of trapped dipolar BECs in the 
TF limit. We explore the static solutions and the low- 
lying collective excitations throughout a large and ex- 
perimentally relevant parameter space, including posi- 
tive and negative dipolar couplings Cdd, positive and 
negative s-wave interactions <?, and cylindrically and 
non-cylindrically symmetric systems. Moreover, our ap- 
proach enables us to unambiguously identify the modes 
responsible for global collapse of the condensate. We 
would like to point out that there is a freely available 
MATLAB implementation of the calculations presented 
in this paper, complete with a graphical user interface, 
which can be found in the supporting material [74]. 

Section [TT| is devoted to the static solutions of the sys- 
tem. Beginning with the underlying Gross-Pitaevskii 
theory for the condensate mean-field, we make the TF ap- 
proximation and outline the methodology for deriving the 
TF static solutions. We then use it to map out the static 
solutions with cylindrical symmetry, for both repulsive 
and attractive s-wave interactions, and then present an 
example case of the static solutions in a non-cylindrically- 
symmetric geometry. We compare to recent experimental 
observations where possible. 

In Section [TIT] we present our methodology for deriving 
the excitation frequencies of a dipolar BEC. This is an 
adaption of the method that Sinha and Castin applied 
to standard s-wave condensates [41] where one consid- 
ers perturbations around the static solutions (derived in 
Section |n| and employs linearized equations of motion 
for these perturbations. At the heart of our approach is 
the exact calculation of the dipolar potential of a hetero- 
geneous ellipsoidal BEC, performed by employing results 
from gravitational potential theory known in astrophysics 
[53H58] and detailed in appendices [B] and [C] 

In Section [TV] we apply this method to calculate the fre- 
quencies of the important low-lying modes of the system, 
namely the monopole, dipole, quadrupole and scissors 
modes, for a cylindrically-symmetric condensate. Wc 
show how these frequencies vary with the key parameters 
of the system, £dd and trap ratio 7, and give physical ex- 
planations for our observations. In Section [V] we extend 
our analysis to non-cylindrically-symmetric BECs. Al- 
though the parameter space of such systems is very large, 
we present pertinent examples. An important feature of 
non-cylindrically-symmetric systems is that they support 
a family of scissors modes which can be employed as a 
test for superfluidity. As such, in Section [VI] we focus on 



these scissors modes and show how they vary with key 
parameters. Finally, in Section VII we summarise our 
findings. 

There are three appendices included in this paper. Ap- 
pendix [A] contains a plot of the frequencies of the collec- 
tive modes of the BEC as a function of Sdd- Appendix [B] 
outlines the method by which we calculate dipolar poten- 
tials due to arbitrary polynomial density distributions of 
atoms. This is the main technical advance of this work 
over our previous papers which were limited to the dipo- 
lar potentials associated with strictly paraboloidal den- 
sity distributions, i.e. those of the same symmetry class 
as the static solution. In Appendix [C] we give a closed 
formula in terms of elliptic integrals for the dipolar po- 
tential inside a triaxial ellipsoid with a parabolic density 
profile. This is a special but important case of the general 
theory outlined in Appendix [B| 



II. STATIC SOLUTIONS 

A. Methodology for obtaining static solutions 

At zero temperature the condensate is well-described 
by a mean-field order parameter, or "wave function" , 
if)(r, t). This defines an atomic density distribution via 
n(r,t) = |^(r,t)| 2 . Static solutions, denoted by ip(r), 
satisfy the time-independent Gross-Pitaevskii equation 
(CPE) given by 0, 



2m 



V 2 + ^(r) + $ dd (r)+. 9 $(r) 



4>(r) = Mr)(4) 



where \i is the chemical potential of the system. The 
external potential V(r) is typically harmonic with the 
general form, 

V(r) = X - m uj\ [(1 - e)x 2 + (1 + e)y 2 + 7 V] . (5) 

Here wj_ is the average trap frequency in the x — y plane 
and the trap aspect ratio 7 = lo z /lj±_ defines the trapping 
in the axial (z) direction. The trap ellipticity in the x — y 
plane e defines the transverse trap frequencies via uj x = 
\/l — e u>± and uj y = \/l + e 10 ± . When e = the trap is 
cylindrically symmetric. 

The $ d d-term in Eq. Q is the mean-field potential 
arising from the dipolar interactions 



$dd(r) 



ra(r')L/ dd (r-r')dV 



(6) 



This term is a non-local functional of the density and is 
the source of the difficulties associated with theoretical 
treatments of dipolar BECs: it turns the GPE into an 
intcgro-differential equation. A key feature of the ap- 
proach taken by us in this paper is to calculate this term 
analytically. To this end we express the dipolar mean- 
field in terms of a fictitious electrostatic potential (j)(r) 
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mum geo 



*dd(r) - ~C dd ( ^0(r) + -n(r) ) , (7) 



where 



are given in Appendix [C] Note that for a cylindrically- 
symmetric trap e = 0, the static condensate profile is also 
cylindrically-symmetric with aspect ratio k x = n y =: k. 
In the cylindrically-symmetric case the integrals fiijk of 
Eq. (12) can be evaluated in terms of the 2-F1 Gauss hy- 



pergeometric function |59J [60J for any i,j,k 



1 r n ( r 'l d 3 r > 



(8) 



ftijk 



2.F1 (k+ l,l;i+j + k+ |;1-k 2 ) 
(1 + 2i + 2j + 2fc)K 2 ( I+ J) 



(13) 



0(r) satisfies Poisson's equation V 2 </i(r) = — n(r). Note 
that in ^ we have taken the dipoles to be aligned along 
the z-direction. The term n(r) /3 appearing on the right 
hand side of ((71) cancels the Dirac delta function which 

. This means that 



arises in the cr(f)(r)dz term [51] 



$dd includes only the long-range (r -3 ) part of the dipolar 
interaction, exactly as written in Equation (pj). 

We assume the TF approximation where the zero-point 
kinetic energy of the atoms in the trap is neglected. 
Dropping the relevant V 2 -term in Eq. Q leads to 



For the parabolic density profile of Eq. (10), the TF 
Eq. ([9| becomes 



/i = 3ge dd 
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+ V(r) + (l-e dd ) 
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R r — 



Pony 
y 2 



3/3oo2^ 2 ] 
(14) 



Inspection of the coefficients of x 2 , y 2 and z 2 leads to 
three self-consistency relations, given by 



V(r) + <l> dd (r)+gn(r) = /i. 



(9) 



For an s-wave BEC under harmonic trapping, the exact 
density profile in the TF approximation is known to be 
an inverted parabola [3] with the general form, 



n(r) = n 1 



V 
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for n(r) > (10) 



R, = 



2gn 
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(15) 
(16) 
(17) 



where n = l5N/(8TrR x R y R z ) is the central density, and Solving Eqs. @-([T7} gives the exact static solutions of 



R x ,R y , and R z are the condensate radii. In order to 
obtain the dipolar potential arising from this density dis- 
tribution, one must find the corresponding electrostatic 
potential of Eq. ^ . References [25l [26] follow this pro- 
cedure, and arrive at the remarkable conclusion that the 
dipolar potential $dd is also parabolic. Therefore, a 
parabolic density profile is also an exact solution of the 
time- independent TF Equation (|9) even in the presence 
of dipolar interactions. In Section In] and Appendices [B] 
and [C] we point out that this result can be extended us- 
ing results from 19th century gravitational potential the- 
ory [531 1541 157j to arbitrary polynomial densities yielding 
polynomial dipolar potentials of the same degree. For 
the parabolic density profile at hand, the internal dipo- 
lar potential is given by [26l E] 

* / n, / x . 3ge dd n K x K y 
$dd(r) = -ge dd n(r) + y - 

x [Ami - {Pioix 2 + Pony 2 + 3/W 2 ) #J 2 ] (11) 

where k x = R x /R z and n v = R y /R z are the aspect ratios 
of the condensate, and 



the system in the TF regime. 

The energetic stability of the condensate is determined 
by the TF energy functional 



E 



V{r) + i$ d d(r) + ^n(r) ) n(r)d 3 r. ( i s ) 



Inserting the parabolic density profile ( 10 ) yields an en- 
ergy landscape 



E = 



15N 2 g 
28Trn x n v R3 



[(l-£dd) 



(7/3ooi " 3A 



002 



K x Pioi 



K 2 /?01l) 



N 



(19) 



ds 



{K 2 X +.S) 1 +I( K 2 +Sp'+5(1 + S )fc+ 



t, (12) 



Static solutions correspond to stationary points in the 
energy landscape. If the stationary point is a local mini- 
mum in the energy landscape, it corresponds to a physi- 
cally stable solution. However, if the stationary point is a 
maximum or a saddle point, the corresponding solution 
will be energetically unstable. The nature of the sta- 
tionary point can be determined by performing a second 
derivative test on Eq. ( 19 1 with respect to the variables 



where k are integers. Explicit expressions for 
/?ooi, /3ioi: /^on, an d /?oo2 in terms of elliptic integrals 



K x ,n y , and R z . This leads to 6 lengthy equations that 
will not be presented here. Note that this only deter- 
mines whether the stationary point is a local minimum 
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within the class of parabolic density profiles. In other 
words, with the three variables K x ,n y , and R z we are 
only able to determine stability against "scaling" fluctu- 
ations, so named because they correspond to a rescaling 
of the static solution [l)TJ ■ However, the class of scal- 
ing fluctuations includes important low-lying shape os- 
cillations such as the monopole and quadrupole modes. 
Although higher order (beyond quadrupole) modes can 
become unstable in certain regimes, as a criterion of sta- 
bility we will use the local minima of ( 19 ). This assump- 



tion is supported by the recent experiments by Koch et 
al. [TT] , where dipolar BECs were produced with e^d > 1 
that were stable over significant time-scales. 



B. Cylindrically-symmetric static solutions for 
g > 0, and the critical trap ratios 7+^ and J~ llt 

We have obtained the static solutions for 



cylindrically-symmetric BEC by solving Eqs. ( 15 ) to ( 17 ) 



numerically. The solutions behave differently depending 
on whether the s-wave interactions are repulsive or at- 
tractive. We begin by considering the g > case. The 
ensuing static solutions, characterised by their aspect ra- 
tio k, are presented in Fig. [2] as a function of £dd with 
each line representing a different trap ratio 7. While the 
TF solutions in the regime £dd > have been discussed 
previously [UJ US] > the regime of £dd < has not been 
studied. Be aware that when we fix g > 0, the regime 
£dd < (left hand side of Fig. [2]) corresponds to Cdd < 
where the dipolar interaction is reversed, repelling along 
z and attracting in the transverse direction. This can be 
achieved by rapid rotation of the field aligning the dipoles 
about the z-axis |16j . 

Before we examine the question of stability, let us first 
interpret the structure of the solutions shown in Fig. [5] 
Imagine an experiment in which the magnitude of £dd is 
slowly increased from zero. At £dd = we have purely 
s-wave interactions and all solutions have the same as- 
pect ratio as the trap, i.e. k = 7. As £dd is increased 
above zero k decreases so that k < 7 for all solutions. 
This is because standard magnetostriction causes dipolar 
BECs to be more cigar-shaped than their s-wave coun- 
terparts. Conversely, if £dd is made negative then k in- 
creases so that k > 7 for all solutions. This is because 
when Cdd < we have non-standard (reversed) magne- 
tostriction which leads to a more pancake shaped BEC. 

Consider now the stability of the solutions, beginning 
with the range — 1/2 < £dd < 1 [white region in Fig. Q]. 
We find that the energy landscape ( 19 1 has only one sta- 



tionary point, namely a global energy minimum, and it 
occurs at finite values of the radii R x (= R y ), and R z . 
This global minimum persists for all trap ratios (out- 
side of the range — 1/2 < £dd < 1 the existence of sta- 
ble static solutions depends on 7). Thus, in the range 
— 1/2 < £dd < 1 the static TF solution is stable against 
scaling fluctuations. Other classes of perturbation could 
lead to instability, but there is good reason to believe 




FIG. 2: (Color online) Aspect ratio k, of the g > 
cylindrically-symmmetric static solutions as a function of Edd 
according to Eqs. ( |15[ )-(|17[). Note that e c id < corresponds 
to Cdd < 0. The solid lines indicate the static solutions for 
specific trap ratios 7 which are equally spaced on a logarith- 
mic scale in the range 7 = [0.1, 10], with black/red lines cor- 
respond to minimum/saddle points in the energy landscape. 
The parameter space of global, metastable and unstable so- 
lutions is denoted by white, light grey and dark grey regions, 
respectively. 



that in this range the parabolic solution is stable against 
these too. Take, for example, phonons, i.e. local density 
perturbations. These have a character that can be con- 
sidered opposite to the global motion involved in scaling 
oscillations. The local character of phonons means that 
considerable insight can be gained from the limiting case 
of a homogeneous dipolar condensate. The energy of a 
plane wave perturbation (phonon) with momentum p is 
given by the Bogoliubov energy Eb [12] . 



El = 



(2 \ 2 2 
jy +2< ? n{l + £ dd (3cos 2 0-l)}|-, (20) 



where is the angle between the momentum of the 
phonon and the polarization direction. The perturba- 
tion evolves as ~ exp(iE&t/K) and so when E^ < the 
perturbations grow exponentially, signifying a dynami- 
cal instability. Dynamical stability requires that E^ > 
which, for g > 0, corresponds to the requirement that 
[1 + £ d d(3cos 2 6» - 1)] > in Eq. pol). This leads once 



again to precisely the stability condition — 1/2 < £dd < 1- 
Outside of the regime — 1/2 < £dd < 1 the global energy 
minimum of the TF system is a collapsed state where at 
least one of the radii is zero, just like in the uniform dipo- 
lar BEC case. However, unlike the uniform case, in the 
presence of a trap the energy functional can also support 
a local energy minimum corresponding to a metastable 
solution [light grey region in Fig. |2J. The existence of 
a metastable solution means there must also be a sad- 
dle point connecting the metastable solution to the col- 
lapsed state and this is indicated by the dark grey region 
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in Fig. @. 

In general, the occurence of metastable solutions de- 
pends sensitively on e^d an d 7- Remarkably, however, 
there are two critical trap ratios, 7* it = 5.17 and 7~ it = 
0.19, beyond which the BEC is stable against scaling fluc- 
tuations even as the strength of the dipolar interactions 
becomes infinite. First consider e^d > 1, f° r which there 
is a susceptibility for collapse towards an infinitely nar- 
row line of end-to-end dipoles (R x = R y — > 0) . Providing 
7 > 7^ it , i.e. if the trap is pancake enough, condensate 
solutions metastable against scaling fluctuations persist 
even as £dd — > oo [T3J HU |2B]. Referring to Fig. [2j these 
curves are located in the upper right hand portion of 
the plot and asymptote to horizontal lines as £dd is in- 
creased (see Fig. 3 in for a plot which extends £dd 
to much higher values than shown here so that this be- 
havior is clearer). However, if the trap is not pancake- 
shaped enough, i.e. 7 < 7„ it , then as £dd is increased 
from zero the local energy minimum eventually disap- 
pears and no stable solutions exist. Referring again to 
Fig. [2j these are the curves that turn over as £dd is 
increased, and in so doing enter the dark grey region. 
Second, consider £dd < —0.5, for which the system is 
susceptible to collapse into an infinitely thin pancake of 
side-by-side dipoles (R z — > 0). If the trap is sufficiently 
cigar-like with 7 < 7~ it collapse via scaling oscillations 
is suppressed even in the limit £dd —00. These curves 
are located in the lower left hand portion of Fig. [2] and 
asymptote to horizontal lines. However, if the trap is 
not cigar-shaped enough, i.e. 7 > 7~ it , then for suffi- 
ciently large and negative £dd the metastable solution 
disappears, bending upwards to enter the dark grey re- 
gion on the left hand portion of Fig. [2] and the system 
becomes unstable to collapse. 

In a recent experiment Lahaye et al. [10] measured the 
aspect ratio of the dipolar condensate over the range 
0^£dd~lj using a Feshbach resonance to tune g, and 
found very good agreement with the TF predictions. 
Similarly, Koch et al. [TJJ observed the threshold for col- 
lapse in a 7 = 1 system to be £dd ~ 11, in excellent 
agreement with the TF prediction of £dd = 1-06. Using 
various trap ratios, it was also found that collapse became 
suppressed in flattened geometries and the critical trap 
ratio was observed to exist in the range 7^ it rs 5 — 10, 
which is in qualitative agreement with the TF predic- 
tions. 



C. Cylindrically-symmetric static solutions for 
g < 0, and the nature of dipolar stabilization 

We now consider the case of attractive s-wave interac- 
tions g < 0. Negative values of g can be achieved using 
a Feshbach resonance, as implemented in a 52 Cr BEC in 
[TT] . The static solutions are presented in Fig. [3] Be 
aware that because g < 0, £dd < (£dd > 0) now cor- 
responds to Cdd > (Cdd < 0). The stability diagram 
differs greatly from the g > case and, in particular, no 
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FIG. 3: (Color online) Aspect ratio k, of the g < 
cylindrically-symmetric static solutions as a function of £dd- 
Note that the regime of £dd > corresponds to Cdd < 0. The 
lines denote static solutions for specific trap ratios 7, equally 
spaced on a logarithmic scale in the ranges 7 = [0.010, 7~ it ] 
(lower right set of curves) and 7 = hi t , 100] (upper left set of 
curves). Arrows indicate direction of increasing 7. The light 
grey region and black lines correspond to minimum points, 
while red lines correspond to saddle points in the energy land- 
scape. At the extreme left and right hand sides of the fig- 
ure the stable solutions become horizontal lines as they tend 
asymptotically to the trap aspect ratio k — > 7 (see text). 



TF solutions exist in the range —1/2 < £dd < 1- Nev- 
ertheless, TF solutions can exist outside of this range in 
regions of parameter space determined by the two critical 
trap ratios 7~ it and 7^ it introduced in the previous sec- 
tion. We find that for £ d d > solutions only exist for sig- 
nificantly cigar-shaped geometries with 7 < 7 CTit = 0.19, 
while for £dd < solutions only exist for significantly 
pancake-shaped geometries with 7 > 7^ it = 5.17. Fur- 
thermore, the attractive s-wave interactions always cause 
the global minimum to be a collapsed state. This means 
that static solutions are only ever metastable (light grey 
region in Fig. [3]) . 

Again, valuable insight can be gained by considering 
the Bogoliubov spectrum (20 1, this time with g < 0. 



Firstly, for the purely s-wave case we recall the well- 
known result [3] that a homogeneous attractive BEC is al- 
ways unstable to collapse. With dipolar interactions the 
uniform system is stable to axial perturbations (9 = 0) 
for £dd < —1/2 and to radial perturbations (8 = tt/2) 
for £dd > 1- This is the exact opposite of the g > case 
and corroborates the lack of solutions given by the TF 
equations for —1/2 < £dd < 1- Of course, £dd < —1/2 
and £dd > 1 cannot be simultaneously satisfied and so 
a uniform dipolar system with g < is always unsta- 
ble. However, when the system is trapped the conden- 
sate can be stabilized even in the TF regime. The mean 
dipolar interaction depends on the condensate shape and 
can become net repulsive in cigar-shaped systems when 
£dd > (for which Cdd < 0), and in pancake-shaped sys- 
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terns when e^d < (for which Cdd > 0). Remarkably, in 
these cases it is the dipolar interactions that stabilize the 
BEC against the attractive s-wave interactions and lead 
to the regions of metastable static solutions observed in 
Fig. [3] Without the dipolar interactions the BEC would 
collapse. 

Although our model predicts that no solutions exist 
for —1/2 < e^d < 1, it is well-known that stable con- 
densates with purely attractive s-wave interactions can 
exist. Zero-point motion of the atoms (ignored in the TF 
model) induced by the trapping potential stabilises the 
condensate up to a critical number of atoms or interac- 
tion magnitude [3]. One can expect, therefore, that for 
a finite number of atoms the presence of zero-point mo- 
tion enhances the stability of the condensate beyond the 
TF solutions. Koch et al. [TT] have produced a dipolar 
condensate with g < and reported the onset of collapse 
for £dd ^ - 7 in a trap with 7 = 10. For this trap the 
TF static solutions disappear for £dd <^ — 1-5. The inclu- 
sion of zero-point motion cannot explain this discrepancy 
between theory and experiment since it should increase 
the critical value of £dd beyond —1.5, not decrease it. 
Furthermore, including the zero-point motion by using a 
gaussian ansatz leads to an almost identical prediction 
[TT] . One possible explanation of the discrepancy is that 
the dominant dipolar interactions may lead to significant 
deviations of the density profile from a single-peaked in- 
verted parabola/gaussian profile, for example, Ronen et 
al. [25] have predicted bi-concave density structures, al- 
beit in the different regime of £dd — > 00. 

The metastable TF solutions shown in Fig. [3] have a 
counter-intuitive dependence upon £dd- Take, for exam- 
ple, the family of metastable solutions (black curves) in 
the lower right hand portion of the figure. We see that 
as £dd increases k decreases (condensate becomes more 
cigar-shaped). This is in contradiction to what one might 
naively expect because on this side of the figure Cdd < 0, 
and so the dipolar interaction has an energetic preference 
for dipoles sitting side-by-side not end-to-end! In order 
to appreciate what is happening in this region of Fig. 
[3j observe that for each value of £dd there is a critical 
value of the condensate aspect ratio k below which the 
system is metastable, and above which it is unstable. As 
£dd is increased from this point the net repulsive dipo- 
lar interactions favor elongating the BEC so that atoms 
sit further from each other, thereby lowering the inter- 
action energy and decreasing k. In the limit n — > one 
can show that the dipolar mean-field potential tends to 
<f>dd = ~.9£dd^( r ) [33]; i.e. it behaves like a spherically- 
symmetric contact interaction which is repulsive when 
g < and £dd > 0. This means that when £dd is in- 
creased in a strongly cigar-shaped configuration the con- 
densate aspect ratio tends asymptotically towards that of 
the trap n — > 7, as it must for a system with net-repulsive 
spherically-symmetric contact interactions. This behav- 
ior can be seen in Fig. [3] where the black curves all tend 
to straight lines as £dd is increased, and the asymptotic 
value of k they tend to is exactly the trap aspect ratio 7. 



A parallel argument holds for the upper left hand por- 
tion of Fig. [3] where the condensate is quite strongly 
pancake-shaped (7 > 7„ it ): in the limit n — > 00 one 
can show that the dipolar mean-field potential tends to 
$dd = 2<7£dd'ra(r) [S3], i-e. it behaves like a spherically- 
symmetric contact interaction which is repulsive when 
g < and £dd < 0. In this portion of the figure one 
therefore also finds that as |£dd| — > 00 the condensate as- 
pect ratio tends asymptotically towards that of the trap 
k — > 7. 

It is tempting to conclude that the collapse that occurs 
as the strength of the dipolar interactions is reduced rel- 
ative to the s-wave interactions is an "s-wave collapse" 
of the type encountered in BECs with attractive purely 
s-wave interactions, which typically occur through an un- 
stable monopole mode |63j . However, from Fig. [3] we see 
that the magnitude of the dipolar interaction is always 
finite at the collapse point. Furthermore, we shall find in 
subsequent sections that it is always a quadrupole mode 
that is responsible for collapse in a TF dipolar BEC. Col- 
lapse via a quadrupole mode has a ID or 2D character, 
depending on the sign of Cdd [31] . and is distinct from 
collapse via the monopole mode which has a 3D charac- 
ter. 

Having indicated how the static solutions behave for 
attractive s-wave interactions g < 0, for the remainder of 
the paper we will concentrate (although not exclusively) 
on the more common case of repulsive s-wave interac- 
tions. 



D. Non-cylindrically-symmetric static solutions 

We now consider the more general case of a non- 
cylindrically-symmetric system for which the trap ellip- 
ticity e is finite and k x and n y typically differ. Note that 
we perform our analysis of non-cylindrically-symmetric 
static solutions for repulsive s-wave interactions g > 0. 
In Fig. [4] we show how k x and k v vary as a function 
of £dd in a non-cylindrically-symmetric trap. Different 
values of trap ratio are considered and generic qualita- 
tive features exist. The splitting of k x and K y is ev- 
ident, with k x shifting upwards and n y shifting down- 
wards in comparison to the cylindrically-symmetric solu- 
tions. Furthermore, the branches become less stable to 
collapse. For example, for 7 = 0.18 < 7~ it (Fig. j4^a) ) , 
in the cylindrically-symmetric system there exist stable 
solutions for £dd — > —00, but in the anisotropic case, sta- 
tionary solutions only exist up to £dd — — 11. 

We already noted in the introduction that for a cylin- 
drically symmetric dipolar BEC magnetostriction causes 
the radial vs axial aspect ratio k = R x /R z to differ from 
the trap ratio 7, in contrast to a pure s-wave BEC for 
which k = 7. It is therefore interesting to note that we 
find that when the trap is not cylindrically symmetric a 
dipolar BEC also has an ellipticity in the xy-plane which 
differs from that of the trap, although the deviation is 
generally small. This occurs despite the fact that dipolar 
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FIG. 4: Stable static solutions, characterised by the aspect 
ratios K x (dotted lines) and k v (dashed lines), in a non- 
cylindrically-symmetric trap with ellipticity e = 0.75 and (a) 
7 = 0.18, (b) 7 = 0.333, (c) 7 = 3 and (d) 7 = 5.5. Stable 
(unstable) static solutions are indicated by black (grey) lines. 
The corresponding static solutions for e = are indicated by 
solid lines. 



interactions are radially symmetric. 



III. CALCULATION OF THE EXCITATION 
SPECTRUM 

Now that we have exhibited some of the features of the 
static solutions in the TF regime, we wish to determine 
their excitation spectrum. The methods which have been 
previously used for finding the excitation spectrum of a 
dipolar BEC include: i) A variational approach applied 
to a gaussian approximation for the BEC density pro- 
file P31 [2H S3 [M] . This allows one to derive equations 
of motion for the widths of the gaussian. ii) Using the 
equations of dissipationless hydrodynamics, namely the 
continuity and Euler equations, to obtain equations of 
motion for the TF radii [3S] [50] . This method is exact 
in the TF limit (recall that the TF regime is mathemat- 
ically identical to the hydrodynamics of superfluids at 
zero temperature), hi) Solving the full Bogoliubov equa- 
tions [5711351131]. iv) Solving for the time evolution of the 
full time-dependent GPE under well-chosen peturbations 

[HI21I37]. 

Methods i) and ii) are simple but yield only the three 
lowest energy collective modes (the monopole and two 
quadrupole modes). However, in the pure s-wave case 
these methods do have the advantage of giving analytic 
expressions for the frequencies, and in the dipolar case 
the frequencies are given by the solution of the algebraic 
equations (15 17), which are simple to solve. This is to be 



contrasted with the other methods which, although more 
general, require much more sophisticated numerical ap- 



proaches. Furthermore, the non-local nature of the dipo- 
lar interactions make numerical calculations considerably 
more intensive than their s-wave equivalents. Therefore, 
the approach we adopt here is semi-analytic, incorporat- 
ing analytic results for the non-local dipolar potential, 
thereby reducing the problem to the solution of (local) 
algebraic equations. 

In our approach we generalize the methodology previ- 
ously applied by Sinha and Castin [3T] to pure s-wave 
BECs, where linearized equations of motion are derived 
for small perturbations about the mean-field stationary 
solution. One strength of this method, in contrast to 
some of those mentioned above, is that it is trivially 
extended to arbitrary modes of excitation and unstable 
modes/dynamical instability. For example, extension of 
the variational approach to higher-order modes (e.g., to 
consider the scissors modes of an s-wave BEC [66]) re- 
quires that this is "built-in" to the variational ansatz 
itself. We outline our approach below. 

The dynamics of the condensate wave function ip(r, t) 
is described by the time-dependent Gross-Pitaevskii 
equation, 



2m 



V + $ 



dd 



V>, (21) 



where, for convenience, we have dropped the arguments 
r and t. By expressing ip in terms of its density n and 
phase S as, 

tp = Vne iS , 



one obtains from Eq. (21 1 the well-known hydrodynamic 
equations, 

dn h , 
dt = 



-V ■ (nVS) 



dS 



2m 



|VSf - V-gn- $ 



dd ■ 



(22) 
(23) 



We have dropped the term (?i 2 /2m\/S)V 2 v^ arising 
from density gradients - this is synonymous with mak- 
ing the TF approximation [3]. Note that static solu- 
tions satisfy the equilibrium conditions dn/dt — and 
dS/dt = -fx/h. 

We now consider small perturbations of the density 
and phase, Sn and SS, to static solutions, and linearize 
the hydrodynamic equations (22 23 1. The dynamics of 



' SS' 




' SS' 




= C 




Sn 




Sn 



the perturbations are then described as 

d_ 

dt 
where 

g(l+e dd K)/m 
V • n V 
and the operator K is defined as 

d 2 f <5n(r')dV 



C = 



(24) 



(25) 



(KSn)(r) 



dz 2 



4-7rlr 



Sn(r). (26) 
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The integral in the above expression is carried out over 
the domain where the unperturbed density of Eq. ( 10 ) 



satisfies n > 0, that is, the general ellipsoidal domain 
with radii R x , R y , R z . Extending the integration domain 
to the region where n + Sn > would only add 0(5n 2 ) 
effects, since it is exactly in this extended domain that 
n = 0(8n), whereas the size of the extension is also pro- 
portional to Sn. Clearly, to first order in Sn, the quan- 
tity EddKSn is the dipolar potential associated with the 
density distribtution Sn. To obtain the global shape ex- 
citations of the BEC one has to find the eigenfunctions 
Sn, 5S and eigenvalues A of operator C of Eq. ( [25] ) . For 
such eigenfunctions equation (24) trivially yields an ex- 
ponential time evolution of the form ~ exp(At). When 
the associated eigenvalue A is imaginary, the eigenfunc- 
tion corresponds to a time-dependent oscillation of the 
BEC. However, when A posesses a positive real part, 
the eigenfunction represents an unstable excitation which 
grows exponentially. Such dynamical instabilities are an 
important consideration, for example in rotating conden- 
sates where they initiate vortex lattice formation [4l]l42|. 
However, in the current study we will focus on stable ex- 
citations of non-rotating systems. 

To find such eigenfunctions and eigenvalues we con- 
sider a polynomial ansatz for the perturbations in the 
coordinates x, y, and z, of a total degree v |41j . that is, 



Sn 



p,q,r 



y 



1 y r 



5S 



b pqr x p y"z r , (27) 



where 



v = max {p + q + r} . 

a pqr =£0 



(28) 



All operators in Eq. ( 25 ) , acting on such polynomials of 



degree v, result again in polynomials of the same order. 
For the operator K this property might not be obvious, 
but a remarkable result known from 19 t,l -century gravi- 
tational potential theory states that the integral in Eq. 
(26) evaluated for a polynomial density Sn, yields an- 



other polynomial in x, y, and z. Its coefficients are given 
in terms of the integrals /3ijk defined in Eq. ( 12 ), and the 
exact expressions are presented in Appendix |B| The de- 
gree of the resulting polynomial is v + 2, and taking the 
derivative with respect to z twice yields another poly- 
nomial of degree v again. Thus, operator (25) can be 



rewritten as a matrix mapping between scalar vectors of 
polynomial coefficients. Numerically finding the eigen- 
values and eigenvectors of such a system is a simple task, 
which computational packages can typically perform. 

We present only the lowest-lying shape oscillations cor- 
responding to polynomial phase and density perturba- 
tions of degree v — 1 and v = 2. These form the 
monopole, dipole, quadrupole and scissors modes. These 
excitations are illustrated schematically in Fig. [I] and de- 
scribed below, where we state only the form of the den- 
sity perturbation Sn, since it can be shown that the cor- 
responding phase perturbation SS always contains the 



same monomial terms. Note that a, b, c and d axe real 
positive coefficients. 

• Dipole modes D x , D y and D z : A centre-of-mass 
motion along each trap axis [68J. The D x mode, for 
instance, is characterised by Sn — ±ax. 

• Monopole mode M: An in-phase oscillation of 
all radii with the form Sn — ±a±(bx 2 + cy 2 + dz 2 ). 

• Quadrupole modes Q* v , Qf z and Q\ z : The Qi 
modes feature two radii oscillating in-phase with 
each other (denoted in superscripts) and out-of- 
phase with the remaining radius. For example, the 
Q* y mode is characterised by Sn — ±a ± (bx 2 + 



cy 



dz 2 ). 



• Quadrupole mode Q^: This 2D mode is sup- 
ported only in a plane where the trapping has circu- 
lar symmetry. For example, in the transverse plane 
of a cylindrically-symmetric system the transverse 
radii oscillate out-of-phase with each other, with no 
motion in z, according to Sn = ±a(x ± iy) 2 ■ 

• Scissors modes Sc xy , Sc yz and Sc xz : Shape pre- 
serving rotation of the BEC over a small angle in 
the xy, xz and yz plane, respectively. The Sc xy 
mode is characterised by Sn = ±axy. Note that a 
scissors mode in a given plane requires that the con- 
densate asymmetry in that plane is non-zero other- 
wise no cross-terms exist. Furthermore, the ampli- 
tude of the cross-terms should remain smaller than 
the condensate/trap asymmetry otherwise the scis- 
sors mode turns into a quadrupole mode |37| . 

Note that, in order to confirm the dynamical stability 
of the solution, one must also check that positive eigen- 
values do not exist. We have performed this throughout 
this paper and consistently observe that when Im(A) ^ 
that Re(A) = and that when Im(A) = that Re(A) ^ 0. 
It is also possible to determine excitation frequencies of 
higher order excitations of the BEC by including higher 
order monomial terms. Such modes, for example, play 
an important role in the dynamical instability of rotat- 
ing systems [34] [41] . 

We would like to remind the reader that they can 
download the MATLAB program [74] used to perform 
the calculations described in this section. It includes an 
easy to use graphical user interface. 



IV. EXCITATIONS IN A 
CYLINDRICALLY-SYMMETRIC SYSTEM 

In this section we present the oscillation frequencies of 
the lowest lying stable excitations of a dipolar conden- 
sate in a cylindrically-symmetric trap. Through specific 
examples we indicate how they behave with the key ex- 
perimental parameters, namely the dipolar interaction 
strength £dd and trap ratio 7. Note that we will discuss 
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the scissors modes in more detail in Section IVIl Here we 
will just point out that two scissors modes exist, corre- 
sponding to Sc xz and Sc yz , while the Sc xy mode is non- 
existant due to the cylindrical symmetry of the system. 



A. Variation with dipolar interactions £dd for g > 

In Fig. [5] we show how the collective mode frequen- 
cies vary with the dipolar interactions for the case of 
g > 0. Although it would seem experimentally relevant 
to present these frequencies as a function of £dd, we plot 
them as a function of the aspect ratio k instead. We do 
this for the following two reasons: i) plotting the frequen- 
cies as a function of £dd is problematic since two static 
solutions (metastable local minima and unstable saddle 
points) can exist for a given value of e^d, h) in the crit- 
ical region of collapse at the turning point from stable 
to unstable, the excitation frequencies vary rapidly as a 
function of £ddi but much more smoothly as a function 
of k, and so it is easier to view the behavior as a function 
of k. For completeness we have included the correspond- 
ing plot of the frequencies, but as a function of £dd, in 
Appendix [XJ Also, analytic expressions for the frequen- 
cies of the M and Q\ modes in a cylindrically symmetric 
dipolar BEC in the TF regime can be found in [55] , 

It is worth pointing out that the condensate shape ac- 
counts for a significant part of the physics of these sys- 
tems, and so k is a good variable to work with. For exam- 
ple, in the problem of a rotating dipolar BEC, the critical 
rotation frequency at which a vortex becomes energeti- 
cally favorable is exactly the same as that in a purely 
s-wave BEC providing one corrects for the change in the 
aspect ratio due to the dipolar interactions |67j . However, 
k alone does not contain all the physics. In the case of 
the calculation of the excitation frequencies this is clear 



from Eq. (25) which depends upon both (V-noV^S 1 and 
EddKSn. The former term has a direct dependence upon 
k via the equilibrium density profile na(r), whereas the 
latter term does not. 

We consider three values of trap ratio 7, which fall into 
three distinct regimes: (1) 7 < 7" it , (2) 7crit < 7 < 7 + it 
and (3) 7+ it < 7. Recall that 7^(7^) is the critical 
value above (below) which there exist stable solutions 
for £dd — > +oo(— 00), see also Fig. [2] In each case the 
aspect ratio of the stable solutions exists over a finite 
range k = [k~, k + ]. We will now discuss each regime in 
turn. 
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FIG. 5: (Color online) Excitation frequencies as a function of 
condensate aspect ratio k for a cylindrically-symmetric trap 
with aspect ratio (a) 7 = 0.18, (b) 7=1 and (c) 7 = 5.5. 
Shown are the results for the modes M (orange, circles), D 
(black, stars), Q\ (red, diamonds), Q2 (purple, squares) and 
Sc xz (= Sc yz ) (green, triangles), (d) Static solutions k for 
7 = 0.18, 1 and 5.5. Vertical dashed lines mark the transition 
from stable to unstable for the static solution, and this coin- 
cides with the point at which one of the frequencies tends to 
zero. Vertical dotted lines mark the point at which the static 
solution ceases to exist altogether. 



!■ 7 < 7crit 

In Fig. [5|a) we present the excitation frequencies for 
7 = 0.18 as a function of k. The corresponding static 
solutions are shown as the left hand curve in Fig. J5jd) 
and confirm that the stable static solutions (solid black 
part of curve) exist only over a range of k — 
with k~ ~ 0.03 and re + « 0.25 indicated by vertical 



lines (dashed and dotted, respectively). For k > k + , 
no static solutions exist and so the excitation frequen- 
cies are not plotted beyond this point [dotted vertical 
line in Fig. [5](a) and most left hand dotted vertical line 
in[5](d)]. For k < kT , the static solution is no longer 
a local energy minimum but becomes instead a saddle 
point/maximum that is unstable to collapse [transition 
marked with dashed, vertical line in Fig. p[a) and most 
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left hand dashed vertical line in|5jd)]. Although this so- 
lution is not stable we can still determine its excitation 
spectrum. Crucially, this will reveal which modes are re- 
sponsible for collapse and which remain stable through- 
out. 

Three dipole modes (stars) exist. Dipole modes, in 
general, are decoupled from the internal dynamics of 
the condensate [3] and are determined by the trap fre- 
quencies uj x ,uj y , and u z . This provides an important 
check on our code. For the cylindrically symmetric case, 
lo x = ujy = cj j_ , and hence only two distinct dipole modes 
are visible. For k < k~ the dipole frequencies remain 
constant, indicating the dynamical stability of this mode. 

In general, the remaining modes vary with the dipo- 
lar interactions. Perhaps the key mode here is the 
quadrupole Q\ mode (diamonds). At the point of col- 
lapse the Qi frequency decreases to zero. This is con- 
nected to the dynamical instability of this mode since 
Re(A) > for k < k~ . The physical interpretation of this 
is that the Q\ mode, which comprises of an anisotropic 
oscillation in which the condensate periodically elongates 
and then flattens, mediates the collapse of the conden- 
sate into an infinitely narrow cigar-shaped BEC. In the 
energy landscape picture, this occurs because the bar- 
rier between the local energy minimum and the collapsed 
Rxy — state disappears for n < k~ . Note that the 
link between collapse and the decrease of the quadrupole 
mode frequency to zero has been made in Ref. [21] • The 
Q2 quadrupole mode (squares) decreases to zero, and be- 
comes dynamically unstable, after one passes into the un- 
stable regime as indicated in Fig.^d). The monopole M 
mode (circles) remains stable for k < kT and increases 
with k above this point. 



aspect ratio k. These observations are specific to the case 
of 7 = 1. 



3- 7>7+ t 

In Fig. [5^c) we plot the excitation frequencies for 
7 = 5.5. For k < K~ , no static solutions exist, and 
for k > k + , no stable solutions exist. Here kT w 3.3 and 
k + w 54 (dotted and dashed vertical lines, respectively, 
in Fig.^c) and (d)). 

Again, the dipole modes are constant, while the re- 
maining modes vary with dipolar interactions. Apart 
from the quadrupole Q\ mode, all modes are stable past 
the point of collapse, including the Q2 quadrupole mode. 
The Q\ mode decreases to zero at the point when the con- 
densate collapses to an infinitely flattened pancake BEC, 
which is again consistent with this mode mediating the 
anisotropic collapse. 



B. Variation with dipolar interactions £dd for g < 



We now consider the analogous case but with g < 0. 



As shown in Section II C stable solutions only exist for 
7 > 7+. it = 5.17 and 7 < 7~ it = 0.19, with no stable 
solutions existing in the range 7~ it < 7 < 7„ it - Hence 
we will only consider the two regimes of (1) 7 < 7~ it and 

(2)7>7 c + ri f 



2 - 7 crit < 7 < 7 cri t 

In Fig. [5^b) we present the excitation frequencies for 
7 = 1 as a function of £dd- Since 7~ it < 7 < 7^ it , 
the solutions exist over a finite range of £dd- In terms 
of k, collapse occurs at both limits of its range, i.e., for 
k < k~ and k > n + , where k~ « 0.3 and k + ~ 2.5 
(dashed vertical lines in Fig.[5]jb) and (d)). 

Since the trap is spherically-symmetric, the dipole 
modes (stars) all have identical frequency, i.e. to±. The 
Qi quadrupole frequency (diamonds) decreases to zero at 
both points of collapse, k~ and k + . In the former case, 
this corresponds to the anisotropic collapse into an in- 
finitely narrow BEC, while in the latter case, collapse oc- 
curs into an infinitely flattened BEC. In the low n regime, 
the Q2 quadrupole mode (squares) becomes unstable just 
past the point of collapse, but shows no instability in the 
opposite limit for k > n + . 

It is interesting to note that the monopole mode (cir- 
cles) shows no dependence on k and therefore the dipo- 
lar interactions, in agreement with |25j . Additionally, we 
find that the aspect ratio of the density perturbation re- 
mains fixed at precisely 1 for all values of the condensate 



!■ 7 < 7crit 

In Fig. |6ja) we present the excitation frequencies in 
a highly elongated trap 7 = 0.18. Stable static solu- 
tions exist only for k~ < k < n + where kT 0.25 and 
k + ss 0.29. In this regime we find that all collective fre- 
quencies are purely imaginary and finite, and therefore 
stable. At the critical point for collapse k ss 0.29 the Qi 
mode frequency passes through zero and becomes purely 
real, signifying its dynamical instability. This shows that, 
as for g > 0, the Qi mode mediates collapse and there- 
fore collapse proceeds in a highly anisotropic manner due 
to the anisotropic character of the dipolar interactions. 
The remaining modes do not become dynamically un- 
stable past the critical point, and only vary weakly over 
the range of k shown. It should also be remarked that 
higher order modes with polynomial degree v > 2 also 
become unstable within the range n~ < k < n + where 
no stable parabolic solutions lie, further highlighting the 
metastability of the g < states and confirming the rel- 
evance of the predictions made by the uniform-density 



Bogoliubov spectrum ( 20 ) for a system in the TF regime. 
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FIG. 6: (Color online) Excitation frequencies as a function of 
condensate aspect ratio k for a g < cylindrically-symmetric 
trap with aspect ratio (a) 7 = 0.18 and (b) 7 = 5.5, with cor- 
responding static solutions shown in figure (c). Included are 
the results for the modes M (orange, circles), D (black, stars), 
Qi (red, diamonds), Q2 (purple, squares) and Sc (green, tri- 
angles). Dashed vertical lines indicate the critical point at 
which the stable static solutions turn into unstable ones, dot- 
ted vertical lines indicate endpoints of branches where static 
solutions cease to exist. 



FIG. 7: (Color online) Excitation frequencies in a 
cylindrically-symmetric trap as a function of the trap aspect 
ratio 7 for (a) e dd = 0, (b) e dd = —0.75 and (c) e dd = 1.5. 
Shown are the results for the modes M (orange, circles), D 
(black, stars), Q 1 (red, diamonds), Q2 (purple, squares) and 
Sc (green, triangles). In figures (b) and (c) the frequencies 
for e dd = are included as dashed, gray lines. 



2. 7> 7c l 



Figure |6|b) shows the mode frequencies in a highly 
flattened trap 7 = 5.5, for which stable static solutions 
exist only in the regime kT < k < n + where n~ w 2.7 and 
k + rj 3.3. Similarly, at the point of collapse k ~ 2.7 the 
Qi mode has zero frequency and is dynamically unstable. 
Well below the critical point the Q2 mode frequency also 
becomes zero and dynamically unstable. 



C. Variation with trap ratio 7 



Having illustrated in the previous section how the ex- 
citation frequencies behave for g < 0, from now on we 
will limit ourselves to the case of g > 0. In Fig. [7] we plot 
the excitation frequencies as a function of 7 for various 
values of £dd- A common feature is that the dipole fre- 
quencies scale with their corresponding trap frequencies, 
such that ujd x = ^d v = and ujd z = 7^. We now 
consider the three regimes of zero, negative and positive 
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1. £dd = 

For £dd = stable solutions exist for all 7 and the 
corresponding mode frequencies are plotted in Fig. (7^a) . 
Our results agree with previous studies of non-dipolar 
BECs where analytic expressions for the mode frequen- 
cies can be obtained, see, e.g., [2] and [3]. The Q2 
quadrupole mode has fixed frequency ujq 2 = \/2u>± . The 
scissors mode frequency corresponds to ws Cj! = ws Cs2 = 

+ 7 2 ^_l, and the remaining modes obey the equation 

u 2 = u 2 x [2 + | 7 2 ± X - Vl6 - I672 + 9 7 M , (29) 

where the "+" and "-" solutions correspond to ujm and 
ujq 1 , respectively. 

2. £ d d < 

For e dd = —0.75 (Fig. ^o)) stable solutions, and col- 
lective modes, exist up to a critical trap ratio 7 max k, 
0.56. Beyond that the attractive nature of side-by-side 
dipoles (recall Cdd < 0) makes the system unstable to 
collapse. 

For all of the modes except the Q\ quadrupole mode 
we see the same qualitative behaviour as for the non- 
dipolar case (grey lines) with the modes extending right 
up to the point of collapse with no qualitative distinction 
from the non-dipolar case. The Q\ quadrupole mode, on 
the other hand, initially increases with 7, like the non- 
dipolar case, but as it approaches the point of collapse, 
it rapidly decreases towards zero. Above y nax ; the Q\ 
mode is dynamically unstable. 

3. £ d d > 

For £dd = 1-5 (Fig. ^c)) stable solutions exist only 
above a lower critical trap ratio 7 min sa 2.3. For 7 < 
7 mm the attraction of the end-to-end dipoles becomes 
dominant and induces collapse. Indeed, we find that the 
frequency of the Qi mode passes through zero and is 
dynamically unstable for 7 < 7 min . Above this, the Q\ 
and Q2 frequencies increase towards the limiting values of 
the non-dipolar frequencies of 1.82w_i_ and V2ujj_ because 
in a very pancake-shaped trap the atoms cannot sample 
the anisotropy of the interactions. The remaining modes 
behave qualitatively like the non-dipolar modes for 7 > 

V. NON-CYLINDRIC ALLY-SYMMETRIC 
SYSTEMS AND RELEVANCE TO 
ROTATING-TRAP SYSTEMS 

In this section we will apply our approach to the 
most general case of non-cylindrically-symmetric sys- 



tems. An important experimental scenario where this 
occurs is when condensates are rotated in elliptical har- 
monic traps. This has provided a robust method for gen- 
crating vortices and vortex lattices in condensates (see 
Ref. [35] for a review). Whilst the trap ellipticity in the 
x — y plane is typically small (in most experiments it 
is of the order of a few percent), the rotation accentu- 
ates the ellipticity induced in the condensate. Indeed, 
one can derive effective harmonic trap frequencies for the 
condensate which show that the effective ellipticity can 
be orders of magnitude greater than the static ellipticity 

The Q2 mode can be pictured as a surface wave trav- 
eling around the edge of the condensate. It has a sim- 
ilar shape to the rotating elliptical deformation of the 
trap, and when the trap is rotated at frequencies close 
to that of the Q2 mode then even a perturbatively small 
trap deformation strongly couples to this mode. When 
viewed from the frame of reference rotating with the trap, 
the excitation of the Q2 mode appears as a bifurcation 
of the stationary condensate into a new stationary state 
which mixes in some of the Q2 mode and the condensate 
therefore develops an elliptical shape in the x — y plane. 
For some ranges of rotation speeds this new stationary 
state is in turn dynamically unstable to the excitation of 
higher order modes [3H [35] H5 • This dynamical insta- 
bility disrupts the condensate and is the first step in the 
process by which vortices enter. Although this process is 
complex, the dynamical instability that initiates it is ac- 
curately described within the TF approximation because 
the modes which are initially excited are of sufficiently 
long wavelength. The predictions obtained within the 
TF approximation are in excellent agreement with both 
experiments [H] and numerical simulations of the GPE 
[42l 145] . Although we will not specifically consider rota- 
tion further here, our methodology can be easily extended 
to this scenario [31]. 

As in Section [II D[ we consider finite trap ellipticity e 
in the x — y plane. In Fig. [8] we present the mode frequen- 
cies as a function of ellipticity e for three different exam- 
ples. There are some important generic differences to the 
cylindrical case. Due to the complete anisotropy of the 
trapping potential the dipole mode frequencies (stars) all 
differ, and are equal to the corresponding trap frequen- 
cies uj x = \/l — ewj_, ojy = + ecu± and us z — 7Wj_. 
The monopole mode is present (circles) and its frequency 
increases with e. Strictly speaking the Q2 mode is no 
longer present due to the breakdown of cylindrical sym- 
metry. Instead we find a new Qi mode appearing (up- 
per diamonds) which corresponds to the Q\ z mode for 
£dd > and the Q\ z mode for £dd < 0. The usual 
quadrupole mode Q 3 ^ is also present. The reader is re- 
minded that the superscript in the Qi mode notation 
refers to the in-phase radii, the remaining radius oscil- 
lates out of phases with the other two. Although there 
are actually three permutations of Q%, only two appear 
for any given value of £dd since linear combinations of 
these and the monopole mode can form the remaining 
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FIG. 8: (Color online) Excitation frequencies in a non- 
cylindrically-symmetric trap as a function of the trap ellip- 
ticity e in the x — y plane for the cases of (a) £dd = and 
7=1, (b) £dd = —0.6 and 7 = 0.8 and (c) Edd = 1-25 and 
7 = 2. Shown are the modes D (black, stars), M (orange, 
circles), Q\ (red, diamonds), Sc xy (blue, triangles pointing 
down), Scyz (green, triangles pointing up, upper branch) and 
Sc xz (green, triangles pointing up, lower branch). 



Qi mode. 

We will now consider the specific features for the cases 
presented in Fig. [8] For edd = and 7 = 1 (Fig. |8ja)), 
the solutions are stable right up to t = 1. At this limit 
the x-direction becomes untrapped and this causes the 
system to become unstable with respect to the dipole D x 



mode, as well the Q* y mode which can now expand freely 
along the x-axis. For Edd = —0.6 and 7 = 0.8 (Fig. |8jb)) 
the solutions become unstable to collapse at e k 0.425. 
Only the lower Q* y mode becomes dynamically unstable 
at this point, indicating that it is the mode responsible for 
collapse, which is towards a pancake shaped system. For 
£dd = 1-25 (Fig. |§Jc)) the solutions become unstable to 
collapse at e m 0.45. We again observe that the same Q\ 
mode mediates the collapse, only this time the collapse is 
towards a cigar shaped system. The other modes remain 
stable. 



VI. SCISSORS MODES 

A fundamental question concerning ultracold dipo- 
lar Bose gases is the nature of their quantum state in 
situations when the attractive portion of their interac- 
tions becomes important, such as in cigar-shaped sys- 
tems aligned along the external polarizing field. Some 
time ago [701 IZD it was noticed that, due to exchange 
effects, repulsive interactions favor simple Bosc-Einstein 
condensation (macroscopic occupation of a single quan- 
tum state) over fragmented Bose-Einstein condensation 
(macroscopic occupation of two or more quantum states) 
|72) . The converse is true in the presence of attractive 
interactions. Fragmentation can therefore be potentially 
studied in attractive s-wave condensates which are stabi- 
lized by their zero-point energy. However, the consensus 
seems to be that in those systems mechanical collapse 
of the BEC occurs before significant fragmentation |73j . 
Dipolar interactions, on the other hand, are partially at- 
tractive and partially repulsive, the net balance being 
tunable via the shape of the atomic cloud. A compre- 
hensive investigation of fragmentation in dipolar BECs 
is beyond the scope of the current paper, but below we 
take a step in this direction by calculating the properties 
of scissors modes of a dipolar BEC. 

Experimentally, one of the simplest indicators of 
whether or not an atomic cloud is Bose condensed is to 
examine the momentum distribution following free ex- 
pansion after the trap is turned off [1 . For example, 
according to the equipartition theorem, a gas at thermal 
equilibrium will expand isotropically even if the trap was 
anisotropic. This is not true for a BEC which, due to 
its zero-point energy, expands most rapidly in the direc- 
tion which was most tightly confined. However, for dipo- 
lar BECs the situation is complicated by the anisotropy 
of the long-range interactions which continue to act at 
some level even as the gas expands [51 [S]. Quantized 
vortices are another "smoking gun" indicating the pres- 
ence of a BEC, but these are not easy to controllably 
generate in the cigar-shaped systems which would be 
of primary interest (although they might be useful in 
cases where Cdd < 0, for which pancake-shaped BECs 
have dominant attractive interactions). Furthermore, in 
cigar-shaped systems with Cdd > 0, the rotation speed at 
which a vortex becomes energetically favorable diverges 
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as £dd increases [57]. Scissors modes, on the other hand, 
offer an alternative vehicle for the investigation of super- 
fluidity in dipolar systems which does not suffer from the 
difficulties mentioned above. 

A detailed account of the scissors mode in a pure s- 
wave BEC can be found in [37]. The scissors mode of 
a trapped atomic cloud (thermal or Bose condensed) is 
excited by suddenly rotating the anisotropic trapping po- 
tential over a small angle. Consequently, the atomic 
cloud will experience a restoring force exerted by the 
trap, and provided the angle of rotation is small, it will 
exhibit a shape preserving oscillation around the new 
equilibrium position. The exact response of the atomic 
cloud to the torque of the rotated trapping potential de- 
pends strongly on the moment of inertia of the cloud. 
Since a superfluid is restricted to irrotational flow, it 
will have a significantly different moment of inertia com- 
pared to a thermal cloud. In particular, when the trap 
anisotropy vanishes the moment of inertia of a superfluid 
also vanishes, whereas in a thermal cloud this is not the 
case. The superfluid scissors mode frequency will con- 
sequently approach a finite value, whereas in a thermal 
cloud it will vanish as the trap anisotropy approaches 
zero |37j . A measurement of the scissors mode frequency 
therefore constitutes a direct test for superfluidity (3"ll37|. 
as has been verified experimentally for non-dipolar BECs 

csana. 



(a) 



(b) 



In the following, we will consider the scissors mode to 
be excited by rotating the trapping potential as well as 
the external aligning field of the dipoles simultaneously 
and abruptly through a small angle, such that the con- 
densate suddenly finds itself in a rotationally displaced 
configuration. Three scissors modes now appear due to 
the three distinct permutations of this mode, namely 
Sc xy (triangles pointing down in Fig. [8]), Sc yz (triangles 
pointing up), and Sc xz (triangles pointing up). Clearly, 
from Fig. [SJ the oscillation frequencies of the scissors 
modes are affected by the dipolar interactions. The ef- 
fect of the dipolar interactions is two-fold. Firstly, since 
the dipolar interactions change the aspect ratio of the 
condensate, both the moment of inertia of the conden- 
sate and the torque from the trapping potential acting 
on it will be altered, which consequently will alter the 
oscillation frequency. Secondly, for the Sc xz , Sc yz modes 
there is an additional force present which is related to 
the relative position of the dipoles. This effect is easi- 
est understood when considering a cigar shaped conden- 
sate. When such a condensate is rotated with respect 
to the aligning field, the dipoles are on average slightly 
more side-by-side than in the equilibrium situation. As 
a result, there will be a dipolar restoring force trying to 
re-align the dipoles, which in turn is expected to affect 
the scissors mode frequencies. Figure [9] schematically il- 
lustrates this process for the Sc xz mode. For a pancake 
shaped condensate the effect is opposite. Since the dipo- 
lar interaction potential is rotationally invariant in the xy 
plane, the dipolar restoring force is absent for the Sc xy 
mode. 



t 
t 




FIG. 9: Schematic illustration of the dipolar restoring force 
for the Sc xz mode. When the condensate is rotated with 
respect to the dipole alignment axis z as in situation (b), 
the dipoles will on average be more side-by-side than in the 
aligned case, situation (a). Since this is an energetically un- 
favourable configuration compared to the aligned case, there 
will consequently be a dipolar restoring force present in (b) 
trying to re-align the condensate, illustrated by arrows. 



Explicit expressions for the scissors frequencies can be 
obtained by performing the procedure outlined in section 
|III| analytically, rather than numerically. We start with 
the frequency uj sxy of the Sc xy mode, in which case we 
only expect an influence of dipolar interactions through 
changes in the geometry, and find 



2wi 



K y + K l 



(30) 



where it should be noted that the quantity in brackets 
is precisely the ellipticity of the condensate. The Sc xy 
frequency does not depend explicitly on the strength of 
the dipolar interactions £dd, but merely on the conden- 
sate ellipticity, which is an indication of the absence of 
a dipolar restoring force as discussed above. The con- 
densate ellipticity turns out to be approximately pro- 
portional to the trap ellipticity, where the constant of 
proportionality is dependent on the dipolar interaction 
strength £dd and axial trapping strength 7. As a result, 
the Sc xy scissors frequencies shown in Fig. [8] are (al- 
most) independent of the trap ellipticity for fixed values 
of £dd and 7. Figure 10 'a) shows the Sc xy frequency as 
a function of the axial trapping strength 7, for various 
dipolar interaction strengths £dd- In the presence of dipo- 
lar interactions, the condensate ellipticity deviates from 
the trap ellipticity e (see Section II D ) , and hence the 



Sc xy frequency also changes when dipolar interactions 
are switched on. In the absence of dipolar interactions 
[dashed line in Fig. [lO^a)], the trap and condensate el- 
lipticity are equal and the Sc xy frequency is independent 
of the condensate size, trap ellipticity, as well as the s- 
wave interaction strength |37j . For very prolate (7^1) 
and very oblate (7 3> 1) traps, the dipolar interactions 
become either mainly attractive or mainly repulsive and 
lose their anisotropic character. The dipolar potential 
becomes contactlike and can be renormalized into the 
s-wave interactions (see Section II C or [33]); which do 



not influence the scissors mode frequency. This effect 
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FIG. 10: (Color online) Scissors frequencies as a function of 
the axial trapping strength 7 for various values of £dd, with 
—0.45 < £dd < 0.9 and increasing in the direction of the 
arrow in steps of 0.15. The dashed line indicates Edd = 0. (a) 
Frequency of Sc xy mode for fixed trap ellipticity of e = 0.1. 
For very prolate (7 <C 1) or very oblate (7 > 1) systems, the 
dipolar interactions renormalize into the s-wave interactions 
and ojaxy returns to the non-dipolar value, (b) Frequency 
u) SX z of the Sc xz mode as a function of 7 for a cylindrically 
symmetric trap, scaled to the non-dipolar frequency u4°L 
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edd 



1 - £dd (1 



cK y /3o02) 



(32) 



Here, the quantity edd appears explicitly and as such the 
frequencies depend directly on the strength of the dipolar 
interactions, an effect we attribute to the dipolar restor- 
ing force. Figure 10 'b) shows the above frequencies for a 
cylindrically symmetric trap and as a function of the ax- 
ial trapping strength 7, for various values of Edd- There 
are two distinct effects to be noted. Firstly, when edd > 
(edd < 0) the scissors frequencies go up (down) for cigar 
shaped systems and down (up) for pancake shaped sys- 
tems. This behaviour is consistent with what one would 
expect in the presence of a dipolar restoring force. Sec- 
ondly, for 7 « 1 and 7 > 1 we see that the u sxz fre- 
quency approaches the non-dipolar value again. For the 
Lo sxy frequency this effect could be explained solely by 
the fact that for such values of 7 the condensate aspect 
ratios return to the non-dipolar values. However, for the 
Sc xz and Sc yz modes we have to account for the appar- 
ent vanishing of the dipolar restoring force as well. To 
see why it plays no part here, we have to analyze the 
expectation values of the quantity R — y/x 2 + y 2 + z 2 . 
For 7 < 1 we have (R) ~ (\z\) —> 00, and for 7 > 1 we 
have (R) ~ (y/ ' x 2 + y 2 ) — > 00. 



Although in both cases 



the torque exerted by the dipolar restoring force is pro- 
portional to (R) and in principle approaches infinity, it 
vanishes relative to the other two quantities contribut- 
ing to the scissors frequencies, namely the moment of 
inertia of the condensate and the torque exer ted by the 
trap, which both scale as (R 2 ) [3S]- In Fig. lCTb) this 
behaviour can be observed for the extremal values of 7, 
where the scissors frequencies approach that of the non- 
dipolar case. 



VII. CONCLUSIONS 



is visible in Fig. 10 a) in the form of the scissors fre- 



quency returning to the non-dipolar value for extremal 
values of 7. Finally, we would like to point out a remark- 
able similarity between the scissors frequencies shown in 
Fig. 10 a), and the trap rotation frequencies at which 



the static solution diagram of a rotating dipolar BEC 
shows a bifurcation point, as investigated in reference 
[55] (see figure 1(b) therein). For all values of 7 and edd 
the scissors frequency is precisely twice the bifurcation 
frequency. Presumably, the underlying connection is the 
fact that the scissors mode Sc xy has the same superfluid 
field as a stationary state of a BEC in a rotating trap. 
However, a deeper investigation into the exact nature of 
the relationship is beyond the scope of this paper. 

Turning our attention to the Sc xz and Sc yz frequen- 
cies, the analytical calculation yields 



^sxz I 1 



1 ^ 1 - Edd (l - 5 K x /t 3// 9 102) 

1 - e d d (l - \k x k v Pqq2) 



1 - 2 



(31) 



In this paper we have performed an investigation into 
the static and dynamic states of trapped dipolar Bose- 
Einstein condensates in the Thomas-Fermi regime. We 
have extended our previous work in this area by ex- 
amining new regimes of dipolar and s-wave interactions 
(namely, positive and negative values of C d d and g) , non- 
cylindrically symmetric traps, and different classes of col- 
lective excitation, including the scissors modes. Our ap- 
proach is based upon the analytic calculation of the non- 
local dipolar mean-field potential inside the condensate 
and allows us to calculate the potential due to an ar- 
bitrary polynomial density profile in an efficient man- 
ner. Using this method, we have examined the stabil- 
ity of static states and collective excitations, including 
the behavior of the collective excitations as a function 
of the trap aspect ratio and ellipticity, and as a func- 
tion of the relative strength of the dipolar and s-wave 
interactions. We consistently find that an instability of 
the Qi quadrupole mode mediates global collapse of a 
dipolar BEC whether g > or g < 0. However, there 
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are two critical trap ratios, 7* it = 5.17 and 7~ it = 0.19, 
beyond which the BEC is stable against scaling fluctua- 
tions (monopole and quadrupole excitations) even as the 
strength of the dipolar interaction overwhelms the s-wave 
one, i.e. when e^d ~> ±oo. In the case of attractive s- 
wave interactions (g < 0), where the dipolar interactions 
can stabilize an otherwise unstable condensate, the mag- 
netostriction seems to act counter-intuitively (see Figure 
[3]), although upon closer examination the behavior can 
be explained by understanding how dipolar interactions 
behave in highly confined geometries. 

We have paid special attention to the scissors modes 
because of their sensitivity to superfluidity, which we 
identify as an issue of particular interest in cigar-shaped 
dipolar condensates due to the possibility of fragmenta- 
tion when the attractive part of the dipolar interaction 
becomes significant. Our expressions for the frequencies 
of the scissors modes include a term due to a restoring 
force which is not present in the pure s-wave case, and 
which we identify as arising due to a long-range dipolar 
re-aligment force. 

A freely available MATLAB implementation of the cal- 
culations outlined in this paper, including a graphical 
user interface, can be obtained online [73] . 
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Appendix A: Collective modes frequencies as a 
function of £dd 



In Section |III| we considered the effect of the dipolar 
interactions on the mode frequencies and plotted this as 
a function of k rather than Edd to remove the problem of 
the static solutions being double- valued. However, since 
£dd is a more obvious experimental parameter, we have 
plotted the corresponding frequency plots of Fig. [5j but 
as a function of £dd in Fig. [TTJ 



Appendix B: Calculating the dipolar potential inside 
a heterogenous ellipsoidal BEC 

In this appendix we will concern ourselves with the 
calculation of integrals of the form 
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FIG. 11: Excitation frequencies as a function of £dd for a 
cylindrically-symmetric trap with aspect ratio (a) 7 = 0.18, 
(b) 7 = 1 and (c) 7 = 5.5. Shown are the results for the 
dipole mode D (black, stars), monopole M (orange, circles), 
quadrupoles Qi (red, diamonds) and Q2 (purple, squares), 
and scissors Sc (green, triangles), (d) Static solutions k for 
7 = 0.18 (bottom curve), 1 (center curve) and 5.5 (top curve). 
Stable solutions are marked with a solid line, unstable solu- 
tions ar marked with a dashed (red) line. Dashed vertical 
lines mark the transition point from stable to unstable. 



where the domain of integration is a general ellipsoid with 
semi-axes R x ,R y , R z , and the point r = (x, y, z) is an in- 
ternal point of the ellipsoid. The square brackets indicate 
a functional dependence. In Eqs. (|8| and (26) we need 



to evaluate this integral in order to obtain the fictitious 
electrostatic potential <fi[pijk]{ r ) arising from a particle 
density of the form 



Pijk = x y J z 



(B2) 



rdx'dj/dz', 



(Bl) 



with k nonnegative integers. By taking linear com- 
binations of the general term we can calculate the 
internal dipolar potential created by an arbitrary density 



distribution because Eq. (Bl ) defines a linear integral op- 
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erator acting upon p(r). 

The physically relevant density distributions naturally 
fall into two classes: 



1. The inverted parabola n(r) given by Eq. ( 10 1 which 
corresponds to the static ground state of the BEC. 



2. The excitations Sn(r) and SS(r) given by Eq. (27) 
which can be written as linear combination of terms 
x l y 3 z k . 

Actually these two classes have some overlap because 
certain low lying excitations (monopole, quadrupole, and 
scissors modes), which would otherwise seem to fall into 
class 2, can be described in terms of the parabolic den- 
sity profile of class 1 but with time-oscillating radii (in 
the case of monopole and quadrupole modes), or time- 
oscillating symmetry axes (in the case of the scissors 
modes). In this appendix we give the general theory 
which works for all density distributions of the form (B2 1. 



In Appendix[C]we specialize to the parabolic density dis- 
tributions of class 1 which is a particular case of the gen- 
eral theory and one is able to present the results in terms 
of well known special functions (the elliptic integrals). 

In the context of calculating gravitational potentials 
in astrophysics one encounters exactly the same integrals 
as here and as such, the problem attracted considerable 
interest even in the 19 th century. Among others, signifi- 
cant contributions to the topic were made by MacLaurin, 
Ivory, Green, Poisson, Cayley, Ferrers, and Dyson. A de- 
tailed historical overview can be found, for instance, in 
references [561 158] . However, for our purposes the most 
important contribution came from N. M. Ferrers, who 
showed that the potential (Bl| associated with the den- 



sity (B2| evaluates exactly to a polynomial in the coor- 
dinates and z |53| . As a matter of general interest 
we will outline the method employed by Ferrers to arrive 
at this remarkable result. 

In his 1877 paper, Ferrers first shows how the internal 
(gravitational) potential of an ellipsoid with semi-axes 
R, 



, R y , and R z , with a density of the form 



P(s) 



1,2,3, 



(B3) 



with 



S = l -R^-W-lV ™ 

±\J X J-Vy 

can be calculated, using integration over so-called ho- 
moeoidal shells, to be 

where 



Q = i- 



R 2 x +(t R 2 y + <J R 2 z +a' 



and 



A = 



^(Rl + a)(Rl+a)(Rl+a). 



(B6) 



Homoeoidal shells are shells situated inside the ellipsoid, 
bounded by equidensity surfaces of the density ( B3 1 . Us- 



ing infinitesimally thin homoeoidal shells, the triple in- 
tegral (Bl) can be reduced to a single integral over the 
variable a. For a detailed account on this integration pro- 
cess, see for instance references [26J [55j [56] , or of course 
the original work by Ferrers [53] . Notably, the right hand 
side of Eq. (B5 1 evaluates to a polynomial in x, y, and z. 

Next, Ferrers noted that whenever p = at the bound- 
ary of the ellipsoid, then differentiation of the potential 
with respect to any of the cartesian coordinates, for ex- 
ample x, yields 



d.T 



<f>\p](x,y,z) = 4> 



dp 
d.r 



(x,y,z), 



(B7) 



which can easily be checked with integration by parts. 
Finally, he noted that any monomial, such as pijk de- 
fined in Eq. ( B2 ) , can be expressed by means of a series 
of differential coefficients of powers s m of the variable s 



defined in (B4|, 



Pijk - 



A. m Pi r dxPdyidz r 

m p+q+r<m 



(B8) 



where the Ampqr are constants and whereby it should 
be noted that the order of differentiation never exceeds 
m. By virtue of the latter observation, we can calculate 
the potential of the above density by repeatedly applying 
the step (B7) in the opposite direction, transferring all 



differential operators appearing in (B8|, from inside the 
integral 4>[pijk] to the outside, since at each step we are 
ensured that the density being integrated over contains 
at least a factor of s, and hence is always equal to on 
the boundary. Thus, we arrive at the following result 



<f>[Pijk](x,y,z) = 



AP+q+r 

™Pi r dxPd y«dz r 

rnpqr 



(x,y,z) 



AP+q+r 

V AW — 

Zv mpv dx p dy q dz r 



in which the </>[s m ] terms are known through Eq. (B5|. 
Recalling that the potentials 4>[s m ] are in fact polyno- 
mial in x, y, and z, and hence also any differential quo- 
tient, we can conclude that the potential 4>[pijk] is also a 
polynomial. Crucial point in the above derivation is the 
observation expressed in equation ( B8 ) , that any mono- 



mial can be written as a series of differential coefficients 
of a function of the homoeoidal shell index variable s, 
which is specific to ellipsoids only. In different geome- 
tries, some monomial densities might yield polynomial 
potentials, but in general this is not the case. 

It remains to determine the precise coefficients of this 
polynomial, a task undertaken by F.W. Dyson who found 
a compact and elegant general expression [54 . Through 
the efforts of Ferrers and Dyson, the triple integral of 
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(Bll which depended on the coordinate r, is reduced to 



a finite number of single integrals appearing in the coef- 
ficients of a polynomial only. 

Although Dyson's formula in principle solves the prob- 
lem, it is not particularly suited to numerical computa- 
tion because it still contains differential operators. We 
therefore employ results from a more recent paper by 
Levin and Muratov 57] , which computes the polynomial 
coefficients of the potential 4>[Pijk] explicitly. Levin and 
Muratov make use of generalised depolarisation factors 
defined as 



(2Z-l)!!(2m-l)!!(2n-l)!!- 



2(;+m+n-l) 



2R 



(B9) 



where m,l,n = 0,1,2..., and Pi mn is defined in Eq. ( 12 1 



Next, we write the exponents of pijk = x l y 3 z k as 

i = 2A + S x ,j = 2/i + Sp, k = 2v + 8 V , 

with X,p,v positive integers such that the 5^, 6 V , 5\ are 
either or 1 for k even or odd (respectively), and 
define 

(T = A + /i + Z/+l. 

In the particular case of calculating the dipolar potential 
we are interested in the second derivative with respect 
to 2, rather than the potential <f>[pijk] itself. Using the 
results of Levin and Muratov then, this quantity is given 
by 



pfi 9 /?* W R k a a pa p-q . , \ 2p+S x ,,2 9 +i„ Jr+5,,-2 

^m^^^^wye e r . !iv9 \ +n r ^ ,fc) ' (bio) 

dz 2 p^Oq^O 7~[ {cr-p-q-ry.(2pSx + l){2qS f , + l){2rS l , + l) 



where 



A V v Q T>2l+S x r, 2m + S n n2n+5 v 
■p(i,j,k) _ \ " \ "* \ "* Jlrn n n x n y It z 



and 



Sir, 



(_2y+ / m+n 

(2/)!(2m)!(2n)!' 



Appendix C: Calculating the dipolar potential inside 
an inverted paraboloidal ellipsoid 

In this appendix we calculate the dipolar potential in- 
side an ellipsoidally shaped BEC with a density profile 
of the form 



n(r) = n 1 



x 

rJ 



V 



R 2 J 



n s 



(CI) 



where no is the density in the center of the BEC. This 
is a particular case of the more general density profile 
p(r) = s n discussed in Appendix |B| above, see Eq. (B3). 



Comparing with Eq. ( 10 1 in the main part of the text, we 



see that density profile (CI) corresponds exactly to the 
static solutions discussed in Section [Til Note that the re- 
sults in this appendix generalize those found in Appendix 
A of reference for the cylindrically symmetric case to 
a triaxial ellipsoid. 

In the case n(r) = hqs considered in this appendix, it is 
straightforward to see that the integral (B5 1 for the ficti- 
tious electrostatic potential <p(x, y, z) inside the ellipsoid 



can be expressed as 

_ n a R x R y R z J 1 
~~ 2 \4 

d 2 

+2x 2 y 2 



d{R 2 x )d{R 2 y \ 



d(R% 



+ 2x 2 z 2 



d{Rl) 



+ zJ 



d(R 2 



i) 2 



+2y 2 z 2 



d 2 



d{Rl)d{Rl) 

d 2 



d(R 2 y )d(R 2 ) 3 d(Rl) 2 



i 4 d 2 



!i' „ , , + \^ \I A {R X1 R V1 R Z ) (C2) 



3" d{RT) 2 3 d(R 2 ) 2 



where 



Ia 



00 da 
A"' 



(C3) 



with A defined in Eq. (B6|. For a prolate (cigar-shaped) 
condensate, defined as being when R z > R y > R x , we 
find 



2 / 

I A = — F arccos 

VR^M V 



R.i 

R* 



E>2 d2 

H z - Hy 
Rl-R 2 



(C4) 



where F(6\m) is an elliptic integral of the first kind whose 
properties are well known |60j . In the opposite case of 
an oblate (pancake-shaped) condensate, R y > R x > R z , 
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then 
Ia 



R 2 - R 2 



F I arccos 



Ri, 



-ft,, — Itr. 



Rl-Rl 



(C5) 



The cylindrically symmetric case of R x — R y is given in 
Appendix A of reference [35]. Thus, the problem of cal- 
culating the electrostatic potential <p(r) reduces to one of 
finding derivatives of elliptic integrals, both with respect 
to the argument 9 and the parameter m. To evaluate the 
derivatives of I a needed in Eq. (C2) we shall make use 
of the results 



0_ 
06 

d_ 
dm 



F{6\m) 
F(9\m) 



1 



\J 1 — m sin 2 9 



(C6) 



E(9\m) F(9\m) 



2m(l — to) 2m 
sin 29 



4(1 - m)VT~^ 



man 



where E(9\m) is an elliptic integral of the second kind 
[50] . and 



d_ 
09 
_d_ 
dm 



E(9\m) = Vl-m sin 2 9 
E(9\m) = mm)-F{9\m) 



2m 



(C8) 
(C9) 



When the external polarizing field is aligned along the 
z-axis then the mean-field dipolar potential $dd(r) is 
given by Eq. |7| 



$ dd (r) = -C, 



(1(1 



d 2 1 



Taking <f>(r) from Equation ( C2 ) one has 
n Q R x R y R z 



8 2 
dz z 



+Ay 
+4z 2 



d 

i {-djR 2 

d 2 

2 d(R 2 y )d(R 2 z ) 
d 2 



4x" 



<:) 2 



d{R%)d{Rl 



d(R 



2 \2 



Ia(Rx,R V , Rz)- 



(CIO) 



Thus, the dipolar mean-field potential $dd(r) inside the 
inverted parabola density profile (CI I is itself a quadratic 



function of position (x, y, z), as given in Eq. (ITT 



$dd(r) = -ge dd n(r) 



x [pooi - (Pwix 2 + Pony 2 + 3/W 2 ) Rz 2 ] 



where Ki = Ri/R z , and the coefficients j3ijk defined in 



(C7) Eq. (12 1 can be seen to be 



P001 

Pwi 
Pon 



-2R\ 
ARl 
ARl 







d{Rl) 

d 2 



d(R 2 x )(R 2 
d 2 



lA(Rx,Ry,Rz) (Gil) 

I A {R x ,R yi R z ) (C12) 

Ia(Rx, Ry, Rz) (C13) 



d(R 2 y)(R 2 z) 

4 5 d 2 

^R z Q^ R 2y lA ( Rx ' R y Rz ">- ( C14 ) 



Using (C6 C9| to perform the required derivatives, an 



explicit expression for the dipolar mean-field potential 
$dd can be given in terms of elliptic integrals. In the 
prolate case, when R z > R y > R x , then 
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F arccos [k^ 
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Whilst in the oblate case when R y > R x > F z we have 
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